The information presented here is placed in the public domain, and was written by Doug Burke. The notebook used to create this page is available, and questions can be asked using the GitHub issues page or via Twitter: https://twitter.com/doug_burke.

It will probably make a lot-more sense if you have already read the previous entries in this "series".

K corrections

In today's installment I want to look at the K correction term I mentioned in my first post when calculating the luminosity of a source given a flux and a redshift. However, there's a few things I have to get through first - such as how to represent a spectral model, integration of these models, converting between units - before I can get to the actual K correction. As that leaves me with about one third of a notebook to fill up, I repeat the process, but this time tracking quantities using the dimensional-tf package.

Last time the notebook was run

Spectral models

Before I get to calculating the K correction for a model, I play around with creating a model that represents the spectrum of a source, as it turns out I'm going to want this later. Since I work for the Chandra X-ray Center, I am going to working with models of X-ray emitting sources, concentrating on the 0.1 to 10 keV pass band (I'll explain in a second why I am talking about wavelengths in terms of energy, but the short answer is because Astronomers).

Theoretically, you might expect for spectral models of X-ray sources - that is, the amount of "flux" from a source - to be

1. in SI units

2. a function of wavelength

However, via a combination of history, the fact that in X-ray astronomy we generally deal with the particle nature of light rather than its wave-like nature$^\dagger$, and the need to compare the model to data, X-ray Astrophysical models generally calculate the flux of a source in photon/cm$^2$/s as a function of energy. That is, if $S(E)$ is the source model, then the photon flux, $f$, is defined as

$$ f = \int_{\rm{elo}}^{\rm{ehi}} S(E) \rm{d}E $$

where elo and ehi are the energy band we are interested in (that is, the wavelength range we are using and $E = hc/\lambda$ to convert between wavelength and energy). For sources observable by the Chandra X-ray Observatory, this energy range is close to 0.1 to 10 kilo electronVolts (referred to as keV from now on), which corresponds to 12 to 0.12 nanometers. This means that the units of $S(E)$ are going to be photons/cm$^2$/s/keV, so it is a "flux density". Although not needed here, I will later on need a conversion factor from keV to Joules, so let's write it now.

$^\dagger$ The detectors on Chandra count photons, although some data is taken with a transmission grating in the light path, for which wavelength is the "natural unit", rather than energy.

As the comments indicate, this routine uses the dimensional-tf package, but its use is "hidden", in that both the input and output values of enConv are Double values. As a reminder, the *~ operator from the package converts a value into a type which "knows" about its dimensions, and /~ is its inverse, "extracting" the value whilst checking that the dimensions match, and applying any unit conversion (and, due to the precedence of the operators, x *~ a /~ b is equivalent to (x *~ a) /~ b).

A quick check that this is working: I know that 1 keV should be about $1.602 \times 10^{-16}\ \rm{J}$, or $1.602 \times 10^{-9}\ \rm{erg}$:

Whilst I'm playing around with the dimensional-tf package (more can be found at the end of the notebook), I can also check that I got the conversion between energy and wavelength correct, taking the product of the Plank constant and the speed of light, c, from the Wikipedia page:

So, if 1 keV is 1.2 nm then 0.1 keV is 12 nm and 10 keV is 0.12 nm, which is what I claimed. Success!

In the example I used putStrLn to display the string, using show to convert a value into a string and ++ to append strings. If I had just left it to the notebook to display a string then it would have also displayed the quotes - e.g.

"The wavelength corresponding to ..."

In general I let the notebook display values, but occasionally I'll use putStrLn.

After that little diversion, I can get back to thinking about how I want to represent the spectral models by Haskell code. The obvious interface - at least to me - is a function that accepts a single parameter (an energy, in keV) and returns the flux density for that energy (in units of photon/cm$^2$/s/keV). That is, I can define the following type synonym to represent such a function:

Now, on its own this turns out to be rather useless, since there is no way to select the model or change any parameters of a model! To get around this, I'll expect models to be partially applied, so that a model function will accept the model parameters and then end with the energy at which to evaluate.

Hopefully an example will help. A common model - both because of its simplicity, and because there are many physical processes which behave like this - is a power law, written as (at least, by X-ray astronomers)

$$S(E) = S_0 E^{-\gamma}$$

where $\gamma$ is the slope of the power law$^\dagger$ and a normalisation factor $S_0$, which represents the flux density of the source at 1 keV (when $E$ is in units of $keV$). The power law gets steeper (so more signal at lower energies) as $\gamma$ increases.

$^\dagger$ it is often referred to as the photon index by X-ray astronomers

Using this I can create a function of type Model by specifying the two parameters: for example

which can then be evaluated at one or more energies:

Now, I don't have to use this approach of creating a separate function for each model, since I can also use plMdl directly; that is

and get the same result, but there are times when it's useful.

If I had not given a type signature, then the compiler would have inferred the type

plMdl :: Double -> Double -> Double -> Double

(well, actually it would have come up with something a bit-more general, involving the Floating type class, but let's ignore that for now). Since Model just means Double -> Double, this is the same as the signature I gave - namely Double -> Double -> Model - although I do admit that it can be surprising when looking at a type signature that looks like it takes two arguments when it actually takes three.

Since Model is a type synonym - so that its just a textual substitution, just like a #define statement with the C pre-processor - I haven't gained any type safety by making this substitution. Hopefully it has provided some clarity by better documenting the intent of my code.

Now I want to visualize the model, so I create a function that plots a model over an energy range (with a hard-coded step size of 0.01 keV) and sets up the axis labels (I don't plan to explain how this code works, since it is quite complex; the Chart wiki contains examples and documentation, but I think it's perfectly fine at the moment to just think of this as "and now, the magic of plotting happens"):

Without an explicit type signature - such as the one given above as a comment - the inferred type is rather long:

This says that the function takes two values of the same type (x) and a function which converts from that type to another (y) and creates a Renderable (I'm going to skip over the () part). The constraints on the x type are that it is enumerable (because the values are used in the statement [elo, elo+de .. ehi]), fractional (since de is not an integer), and belongs to the mysterious PlotValue type class, which y also has to. This typeclass is used by Chart to indicate that a type can be included in a plot and, fortunately for me, Double values are an instance of it. I should also add that the choice of x and y as names of types here has nothing to do with X and Y axes of the graph, in case you were wondering.

Hmmm, that's not particularly spectacular now, is it? What I really want is to display the data on a log-log axis (that is, use log scales for both axes). One way to do this is to "wrap" the values to be plotted in the LogValue constructor, which tells Chart to use a log scale. That is, in plotModel, instead of

  cs = map (\e -> (e, mdl e)) es

use

  cs = map (\e -> (LogValue e, LogValue (mdl e))) es

To make this configurable, I create a helper function and then two "plot" routines: plotModel for linear axes and plotModelLL for log axes. I could also create linear-log and log-linear variants, but I'll leave that as an exercise for the interested reader!

The plotModel variant uses the id function, which has an interesting type:

:type id
a -> a

Due to parametricity, this function can only return the value it is sent (the compiler does not know anything about the type, so it can't create or change a value, so the only thing it can do is return the input$^\dagger$), which makes it seem rather pointless. However, it is surprisingly useful; in this particular case it is used where a conversion function is used but no conversion is actually needed.

$^\dagger$ well, it could also return an error or undefined value, but as both of these would end up in a run-time error I'll just look slightly embarassed, mumble a bit, and forget I even mentioned it.

The two versions of the plots for this model look like

How about a model that looks a bit-more interesting, such as a gaussian, defined as

$$S(E) = S_0 e^{-0.5 (E-E_0)^2 / \sigma^2}$$

The plot appears cut off at the high-energy end because I chose an upper limit of 10 keV. Increasing this upper limit leads to more of the Gaussian being drawn:

Compound models can also be created; cMdl also has a type of Model and can be used in the same way as the individual components:

Calculating fluxes

In order to calculate the K-correction for a model I need the flux for a model over a set energy range; that is, I want $\int S(E) dE$. The models I have used so far have analytic solutions for this integral, but instead of actually doing any thinking$^\dagger$, I'm just going to use the same integration technique as in previous notebooks and use an approximation (this time with simpson rather than trap in case the models are not well behaved).

$^\dagger$ also known knowadays as 'using Wolfram Alpha'

Using the power law model, this calculates the 0.1 to 10 keV photon flux to be:

As a quick check, in part to verify that the integration package is working but mainly to show off some simple LaTeX, the analytic solution for the power law model is

$$ \begin{align} \int_{e_1}^{e_2} S_0 E^{-\gamma} {\rm d}E && = && \frac{S_0}{1-\gamma} (e_2^{1-\gamma} - e_1^{1-\gamma}) && \gamma \ne 1 \\ && = && S_0 \rm{log}_e(e_2 / e_1) && \gamma = 1 \end{align} $$

This can be written as (neglecting the fact that I should not be using an equality check for floating-point numbers, so this code is not robust):

So the approximate value calculated by iModel is essentially the same as the exact solution (for the accuracy I need). The exact solution could be calculated for the gaussian model, but I'll leave that for the interested reader (the erf or hmatrix-special packages may be of interest, although I'm sure there's other packages that provide an erf function).

Calculating useful fluxes

Actually, what I really need is not the flux in photon/cm$^2$/s but in erg/cm$^2$/s, that is $\int E S(E) \rm {d}E$. First up, I want a simple way to convert an energy from keV to erg, which is just what enConv, created way back at the start of the notebook, does. Using this, I can modify iModel so that it integrates up energy * photon flux density (i.e. energy flux density) instead:

Note that there's nothing in the types that tell me what the units are; that is, the compiler can't tell me if I've made a mistake and used iFModel when I meant iModel, or used energies in eV rather than keV. It's left to the user - along with the documentation - to ensure things are correct. Once I've got to calculating the K correction, I show how the dimensional-tf package can be used to rewrite this code and ensure that the dimensions match. There are also alternative approaches, such as creating types to represent the values of interest, which are less general but may be easier to use. I may play around with this idea in a future notebook.

With all this, I can compare the two fluxes (photon and energy):

I can compare these results to those from the software that X-ray astronomers use. In this case I use Sherpa, which is a Python-based analysis system. The following commands set up a power-law model with $\gamma=2$, normalisation of 1 and then calculate the photon and energy fluxes over the range 0.1 to 10 keV:

sherpa> dataspace1d(start=0.01, stop=11, step=0.01, id=1)
sherpa> set_source(powlaw1d.plmdl)
sherpa> plmdl.gamma = 2
sherpa> plmdl.ampl = 1
sherpa> calc_photon_flux(lo=0.1, hi=10, id=1)
        9.9000999000999279
sherpa> calc_energy_flux(lo=0.1, hi=10, id=1)
        7.3812306496266669e-09 
sherpa> calc_photon_flux(lo=0.1, hi=10, id=1) / 9.899999999978695
        1.0000100909213367
sherpa> calc_energy_flux(lo=0.1, hi=10, id=1) / 7.3833833576195705e-9
        0.99970843881610427

The numbers are very close$^\dagger$, but it isn't quite a fair comparison, since the Sherpa models are designed to be compared to binned data, which means that the energy flux is approximated by calculating the photon flux within a bin (i.e. range of energies) and then multiplied by the mid-point energy of the bin. That is, the energy flux, $f_{\rm E}$, is calculated using

$$ f_{\rm E} \simeq \sum_{i} \frac{e_{\rm{lo};i} + e_{\rm{hi};i}}{2} s(e_{\rm{lo};i}, e_{\rm{hi};i}) $$

where $s(a,b)$ represents the model integrated over the energy range $a$ to $b$, and the bins cover the desired energy range (in the example above they are defined by the dataspace1d call). In early drafts of this notebook I went off and showed that I get similar results using this approximation scheme, but I think I'll leave that for another time, and get back to the K correction.

$^\dagger$ This level of agreement is significantly better than the accuracy we can get for most Astronomical measurements!

K corrections

The K correction is used in the equation to convert a flux ($f_X$) into a luminosity ($L_X$) for a source at a redshift $z$:

$$ L_X = 4 \pi d_L(z)^2 k(z) f_X $$

where $d_L$ is the luminosity distance and $k(z)$ is the k correction. The reason it is needed is that if you have a source far enough away that you need to worry about including General Relativistic effects when calculating $d_L$, then you also need to account for the fact that - due to the expansion of space - the light we observe has been "redshifted". That is, if a source has a feature at a energy $E_0$ then we would observe that feature at an energy of $E_0/(1+z)$ when the source is at a redshift of $z$. Hopefully an image, created with 1000 characters or so, will be better than 1000 words:

So, although we are observing in the 1 to 5 keV range (in this example), the emission we measure is actually from the 3 to 15 keV range (i.e. $1+z=3$), which will not (for virtually all objects we are interested in) have the same flux. In this case this is indicated by the two different pass bands (i.e. 1-5 and 3-15 keV) having different shapes (and, if you summed up the model) fluxes. If this shift is not corrected for then you can not compare luminosities for objects at different redshifts. What we need is a correction for this, which is where the K correction comes in.

Following Appendix B of Jones et al. 1998, ApJ, 495, 100-114 (or the calculating K corrections document I wrote for people analysing Chandra data), the K correction is calculated as the ratio of the flux of the source in the rest and observed frames, namely

$$ k(z) = \frac{\int_{e_1}^{e_2} E S(E) {\rm d}E}{\int_{e_1 (1+z)}^{e_2 (1+z)} E S(E) {\rm d}E} $$

One way of writing this is to send in the observed frame energy band (elo to ehi), a redshift (z), and a model:

As a quick test:

Interestingly enough, the "default" model I have been using - a power-law with $\gamma=2$ - actually has a K correction of unity here. Is this a fluke, or did I pick a special value? Well, for this model the K correction can be calculated analytically by expanding out $S(E) = S_0 E^{-\gamma}$:

$$ \begin{align} k(z) && = && \frac{\int_{e_1}^{e_2} E^{1-\gamma} {\rm d}E}{\int_{e_1 (1+z)}^{e_2 (1+z)} E^{1-\gamma} {\rm d}E} \\ && = && (1+z)^{\gamma-2} \end{align} $$

So, for a power law, the k correction:

• depends on the slope and redshift, but not the energy range
• the K correction will be 1 when $\gamma-2=0$, i.e. $\gamma=2$ as used in myModel

What does it look like for other models?

The value of the K correction depends, in general, on the redshift of the object, the energy band being used, and the model parameters. It need not be a nice monotonic function, as I show here, when looking at a 3 keV gaussian line with a sigma of 1.2 keV for two different energy bands.

So, with this, I can finally calculate luminosities of sources, assuming I have their redshift and flux.

What would it mean to include dimensions and units?

As I mentioned earlier, using a routine like iFModel, which has a signature of

iFModel :: Double -> Double -> (Double -> Double) -> Double

(where I have expanded out the Model synonym) requires the user to make sure that they understand what units the parameters are. There is no way for the compiler to spot when a user accidentally gives an energy in Joules rather than keV, or a model that calculates an energy, rather than flux, density. I started these notebooks by seeing how easy it was to track quantities and units in Physical calculations, so let's try that approach here by repeating the above steps but using types which track dimensions.

I originally tried to do this with the units package, but ended up getting a bit stuck, so I ended up with dimensional-tf instead. Hopefully I can use the following to guide me when I find time to try again with units.

First I start with defining a bunch of types and functions; they aren't always required$^\dagger$, but provide documentation and simplify some code, and I based the naming scheme on the existing dimensional-tf code.

$^\dagger$ I actually developed the code without explicit types, checked what the compiler had inferred to make sure it made sense, and then created these types.

I am going to use the convention that names that end in P are versions of the "unit-less" code that have been converted to deal with photon fluxes or photon flux densisites. So ModelP is the version of

type Model = Double -> Double

The conversion of the power-law model is relatively easy, but there is extra effort required to deal with the units, in particular with operator precedence. As an example, defining

alpha = -gamma P.*~ P.one

results in a type error because - has a lower precedence than P.*~; that is, the above is taken to mean

alpha = - (gamma P.*~ P.one)

and the error looked like

No instance for (Num (P.Quantity P.DOne Double)) arising from a use of syntactic negation
In the expression: - gamma P.*~ P.one
In an equation for ‘alpha’: alpha = - gamma P.*~ P.one

One thing I do like about the following is that it makes it explicit that the normalization refers to the value at 1 keV:

Evaluating the model returns both the numeric value and units:

Since I'm going to use this set of parameters, I can create an instance of the model in the same way I created myModel above:

Values can be "extracted" from the result either explicitly:

or by using routines such as unFdens (this also shows an example of using different units, in this case eV rather than keV):

As examples of the "type safety" that this approach brings, trying to use a wavelength as an argument causes an error:

I have to explicitly convert the length to an energy:

It's also an error to treat the output value incorrectly. As an example, if I try to calculate an energy flux (by trying to treat the answer as if it has units of W/m$^2$/s) then I also get a compile-time error:

A version of the gaussian model can also be created:

and used:

These models cal also be combined, although you have to remember to use the "correct" mathematical operators (in this case, P.+ rather than just +), otherwise you get a type error:

Unfortunately, because of the types, I can't use my plotModel routine to display this model. It's a relatively simple conversion, but this notebook is long enough already, so I'll leave it for another time.

Given that I can create a model with units of photon/cm$^2$/s/keV (photonFluxDensity), how do I integrate this to get photon/cm$^2$/s (photonFlux)? That is, I want a version of

iModel :: Double -> Double -> Model -> Double
iModel elo ehi mdl = result (absolute 1.0e-6 (simpson mdl elo ehi))

Since the integration code only works with Double values, i.e. it would be a type error to send simpson a model with type ModelP, the code will have to manually convert between types with dimensions and Double values:

So, does this work? The photon flux for myModel over the range 0.1 to 1 keV was 9.9, so what do I get now?

Yay - that's only a factor of $10^4$ out, which is quite close for Astronomy!

The reason for the discrepancy is that the "native" units used by iModelP are SI, so the output is in photon/m$^2$/s rather than photon/cm$^2$/s, which is what iModel returns. This explains the difference of $(100)^2$. Converting the output to my units of choice, using photonFlux, returns the expected value of 9.9 (give or take a little bit of floating-point accuracy):

Fluxes

Using iModelP I can calculate the photon flux, that is, the flux in photon/cm$^2$/s, but what I really need is the energy flux, i.e. erg/cm$^2$/s. First stop: more types and helper functions!

Calculating the energy flux should look familiar:

As a reminder, the photon flux was calculated with

iModelP :: Energy -> Energy -> ModelP -> PhotonFlux
iModelP elo ehi mdlP =
  let x1 = unKeV elo
      x2 = unKeV ehi
      wmdl x = mdlP (keV x) P./~ photonFluxDensity
  in iModel x1 x2 wmdl P.*~ photonFlux

Finally I can calculate the "energy flux" and compare with value from Sherpa:

So I agree with Sherpa, since

The K correction for a model can then be calculated using iModelE:

This is very similar to the original code:

kcorr :: Double -> Double -> Double -> Model -> Double
kcorr elo ehi z mdl = 
  let obs  = iFModel (elo * zp1) (ehi * zp1) mdl
      rest = iFModel elo ehi mdl
      zp1  = 1 + z
  in rest / obs

Does this version give the same answer as the original version?